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ABSTRACT 



Context. The sticking of micron sized dust particles due to surface forces in circumstellar disks is the first stage in the 
production of asteroids and planets. The key ingredients that drive this process are the relative velocity between the 
dust particles in this environment and the complex physics of dust aggregate collisions. 

Aims. Here we present the results of a collision model, which is based on laboratory experiments of these aggregates. 
We investigate the maximum aggregate size and mass that can be reached by coagulation in protoplanetary disks. 
Methods. We use the results of laboratory experiments to establish the collision model (Giittler et al. (2009)). The 
collision model is based on some necessary assumptions: we model the aggregates as spheres having compact and 
porous 'phases' and a continuous transition between these two. We apply this collision model to the Monte Carlo 
method of Zsom & Dullemond (2008) and include Brownian motion, radial drift and turbulence as the sources of 
relative velocity between dust particles. 

Results. We model the growth of dust aggregates at 1 AU at the midplane at three different gas densities. We find that the 
evolution of the dust does not follow the previously assumed growth-fragmentation cycles. Catastrophic fragmentation 
hardly occurs in the three disk models. Furthermore we see long lived, quasi-steady states in the distribution function 
of the aggregates due to bouncing. We explore how the mass and the porosity change upon varying the turbulence 
parameter and by varying the critical mass ratio of dust particles. Upon varying the turbulence parameter, the system 
behaves in a non-linear way and the critical mass ratio has a strong effect on the particle sizes and masses. Particles 
reach Stokes numbers of roughly lO"** during the simulations. 

Conclusions. The particle growth is stopped by bouncing rather than fragmentation in these models. The final Stokes 
number of the aggregates is rather insensitive to the variations of the gas density and the strength of turbulence. The 
maximum mass of the particles is limited to « 1 g (chondrule sized particles). Planetesimal formation can proceed via 
the turbulent concentration of these aerodynamically size-sorted chondrule-sized particles. 

Key words, planets and satellites - formation, accretion, accretion disks, methods: numerical 
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In the c ore accretion pa.rad igm of planet formation (jMizund 
(|1980( ): iPollack et al.l (fToOfii ) planets are the outcome of 
an accretion process that starts with micron-size dust 
grains and covers 40 magnitudes in mass. It can be divided 
into three stages. The first stage of the formation of rocky 
planets and the rocky cores of gas giant planets starts 
with the coagulation of dust in the protoplanet ary disks 
surrou n ding many pre-main-se q uence s tars ([Safronov 
1969r) . IWeidenschilling fc Cuzzil (Il993l ). H um &: WurmI 



2003)). The next stage of planet formation is the for- 



mation of protoplanetary cores from the planetesimals. 
The idea is that the kilometer size planetesimals are 



* This paper is dedicated to the memory of our dear friend and 
colleague Frithjof Brauer (14th March 1980 - 19th September 
2009) who developed powerful models of dust coagulation and 
fragmentation, and thereby studied the formation of planetesi- 
mals beyond the meter size barrier in his PhD thesis. Rest in 
peace, Frithjof. 



SO large, that gravity starts to take over and leads to 
the gravitational agglomeration of these bodies to rocky 
planet s. This scenario was studied already by I Safronov! 
(|19690 . and has since been model ed using numeric al meth - 
ods by [W cidcnschilling (19 80l). iNakagawa et all (Il983l) 
'Mi zuno et all (119881). iSchmitt et al.l (|1997|) . IWetheril 
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(19901), iNomura fc Nakagawal (|2006D . iGaraud fc lT 
(200l), iTanaka et al.l ()2005D and several more authors. 
These models solve for the size distribution of dust ag- 
gregates in the disk as a function of time, and investigate 
if, where and how larger dusty bodies form, and how 
long that takes. Finally, in the third stage, gas accretes 
onto these protoplanets forming giant planets or - in the 
absence of gas - gravitational encounters between these 
protoplanets result in a chaotic, giant impact pha s e, unt il 
orbital stabili t y ha s been achieved ( Chamberi ()2001l ): 
iKokubo et all pOOGl ): iThommes et al.l (|2008l )). 



In this study, we focus on the first phase and address 
the fundamental question of how effective dust growth 
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by surface forces really is; that is, how big do particles 
become by simple sticking processes only. It is known that 
initially, for micron size grains, the growth is driven by 
Brownian motion. This typically leads to sl ow collisions 
and f or ms aggregat e s of fractal structure (jKempf et al.l 
(|1999[ ): iBlum et all (|1996l )). In the current picture of 
dust growth, as these aggregates grow, at some point the 
growth will leave the fractal regime, and collisions will 
start to lead to compac tion and breaking of the aggregates 
(|Blum fc WurmI (|2000l )). embedding of small bodies into 
larger aggregates (leading to 'filling up' of these larger 
aggregates an d compaction due to the force of the collision 
(|Ormel et al.l (p007, )). As the size of t he dust aggregat es 
increases, di fferential vertic al settling (jSafronov J 196911V 
radial dr ift (IWhippld (119721')') an d turbulence (IVolk et al.l 
(|1980D : iMizuno et al.l (|1988D : lOrmel fc Cuzzil (120071 )1 
will become important new mechanisms driving relative 
velocities between aggregates. The increasing relative 
velocities caused by these mechanisms will at least partly 
compensate the lower collision probability due to lower 
surface-over-mass ratio of large aggregates. When the 
aggregates grow to sizes of millimeter to meter, however, 
the st icking efficiency drops strongly (e.g. iBlum fc MiinchI 
(|1993| )) and the relative velocities become so large that 
aggregates can fragment (jBlum fc 'WurmI (I2008D . so called 
'fragmentation barrier'). Another hurdle t hat the particles 
have to circumvent is the 'drift barrier' (jWeidenschillind 
(|l977al l). namely that millimeter, centimeter sized particles 
are lost t o the star d ue to radial drift in a short timescale. 
Recently, lOkuzuniil (j2009l ) pointed out the existence of a 
'charge barrier', which possibly halts the particle growth 
already at an early stage of fractal aggregates. Despite 
many years of efforts, it is not known if the coagulation pro- 
cess can overcome these barriers. These barriers have been 
and still are the main open question of the initial stages 
of planet formation: the growth from dust to planetesimals. 

Several mechanisms have been proposed to overcome 
this p roblem, among w hich arc the trapping of dust i n vor - 
tices (iBarge fc Sommeriai (|19950 : iKlahr fc Hennind (jl997l ). 
iLvra et al.l (l2009l )V trapping of decimeter-sized boulders 
in turbulent eddies and the subsequent g ravitational col- 
lapse of swarms of these trapped boulders (jJohansen et al.l 
(|2007l )). the trapping of particles in a pressure bump cause d 
by the evaporation front of water (|Kretke fc LinI (|2007l ): 
iBrauer et al.l (|2008bD 1 and many more scenarios. However, 
the correct modeling of any of these scenarios requires the 
detailed knowledge of the collisional physics, and these 
models have so far relied either on simplified input phy- 
isics or on simplified initial conditions. 

Because of their complexity, collisional evolution models 
have to make simplifying assumption concerning the out- 
come of dust aggregate collisions, for example that collisions 
always result in sticking, or otherwise use simple recipes 
for the collisional outcome. Ideally, one requires to know 
the detailed outcome of every collision. But modelling this 
microphysics within an evolution model is simply unprac- 
tical. There are computer prog rams that model su c h ind i- 
vidual collisio n s in d etail (e.g. iDominik fc Tieleni (|l997l ). 
ISuvama et al.l (|2008D : Geretshauser et al, in prep.), but 
each model collision takes anywhere from hours to weeks 
to run on a computer. They are therefore not practical 
to use at run-time in a model that computes the overall 
time-dependent evolution of the dust size distribution in- 



side protoplanetary disks. Moreover, such collision models 
themselves often depend on poorly known input physics. 

Another approach to obtain the collisional outcome of 
dust aggregates is to model these collision in the laboratory. 
From the many experiments that have as of now been per- 
formed a picture emerges of the outcome of dust aggregate 
collision under a variety of conditions in the protoplanetary 
disk (PPD). In a companion paper (Giittler et al. (2009), 
henceforth Paper I), we have collected data from over 19 
experiments, and constructed a set of formulae that reason- 
ably well describe the outcomes of these collisions in such a 
way that they can be used as input for models that address 
the temporal evolution of the dust size distribution. 

In this paper we will directly rely on the outcome of 
these laboratory experiments for modeling the dust aggre- 
gate size distribution. As described in Paper I, we have 
produced a mapping of all available collision experiments 
regarding silicate-like particles. The velocity range of these 
experiments is also sufficiently wide to cover various disk 
models which roughly correspond to the conditions at 1 AU 
in the PPD. For details regarding the collisional mapping, 
we will refer to paper I, but we will summarize the elements 
of our new collision model in Sect. 3.1. 

We build this collision kernel into a Monte Carlo code 
for modeling the size- and porosity distribution of dust 
in a protoplanetary disk (Zsom fc DuUcmond (200^), 
hereafter ZsDOS). The outcome of our laboratory-driven 
dust coagulation model is hard to a priori predict since the 
key variables involved depend on a non-trivial interplay 
between the collision kernel (Paper I) and the velocity field. 
We can, however, anticipate two scenarios. In the first, 
particle growth will proceed beyond the meter-size barrier, 
all the way to planetesimals. In the second scenario, growth 
will terminate at an intermediate size. In this case further 
growth to planetesimal sizes may proceed through con- 
centration and subsequent gravit ati onal collapse of the se 
particles (jJohansen et al.l (120071 ). ICuzzi et al.l (|2008l )). 
Thus, our model will provide the starting conditions for 
these concentration models. We do emphasize, however, 
that in this work we do not in any way 'optimize' the 
outcome by laboriously scanning all the parameter space 
or treating environments that may be more conducive for 
growth, like nebula pressure bu mps or trapping o f dust in 
vortices (Kretke fc Lin (2007), iLvra et al.l (|2009l ')'). These 
are obvious expansions of our work. But by considering 
the sensitivity of a few key parameters (e.g., gas density, 
and turbulence strength) on the outcome of the growth 
process, we do obtain a picture of where the arrow of coag- 
ulation typically points to in protoplanetary environments: 
pebbles, boulders or planetesimals. 

In this paper we describe the three nebulae models 
used in this work and the sources of relative velocity be- 
tween the aggregates (Sect. |2|), how we build the coagula- 
tion/fragmentation model of Paper I into the Monte Carlo 
code (Sect. [3]), and what these first results look like (Sect. 
E]). We also test the sensitivity of the results with respect to 
variations in gas density, the velocity field, and other key 
model parameters. Section |5| reflect the importance of our 
result in the context of planetesimal formation and provide 
suggestions for future experiments. Finally, Sect. |5] lists our 
main conclusions. 
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2. The nebulae model 

2.1. Disk models 

In this Section we briefly describe the disk models consid- 
ered in this paper. 

The low density model: Resolved millimeter emission maps 
of protoplanetary disks seem to indicat e a sh allow sur- 
face density profile ([Andrews fc Wilhami (|2007l )): Sg(r) oc 
r""-^. Systematic efi^ects of some of their assumptions, such 
as the disk inclinations or the simplified treatment of the 
temperature distri bution, may s ugg est somewhat steeper 
profiles. Therefore, iBrauer et all pOOSaD adopted the fol- 
lowing profile: 



'^^cm2 (au) 



-0.8 



(1) 



Here we assumed that the central star is of solar mass, the 
disk extends from 0.03 AU until 150 AU and that the total 
mass of the disk is 0.01 M©. Assuming that the pressure 
scale-height is Hp = 0.05 x r and the vertical structure is 
gaussian: 



cxp(-zV2Hp'), 



(2) 



the density at 1 AU in the midplane (z — 0) is 2.4 x 10~^^ 
g cm~'^, approximately two orders of magnitude lower than 
the Minimum Mass Solar Nebulae (MMSN) value. 

MMSN model: The Minimu m Mass Solar N e bulae m odel 
(MMSN) was i ntrod uced by IWeidenschillind (|l977bl ) and 
iHavashi et al.l (jl985l ). From the present state of the Solar 
System today, it is possible to obtain a lower limit to 
the mass if the solar nebulae from which the planets were 
formed. The model assumes that the planets were formed 
where they are currently located (no migration included) . It 
also assumes that all the solid material presented in the so- 
lar nebula had been incorporated in the planets. The loss of 
solid material due to radial drift is not taken into account. 
Despite these uncertainties, the MMSN model is frequently 
used as a benchmark. The surface density of the MMSN 
disk is given by: 



Eg(r) = 1700 



cm2 VAU/ 



-1.5 



(3) 



which corresponds to a total disk mass of 0.01 M© con- 
tained between 0.4 and 30 AU (between the orbits of 
Mercury and Neptune). Assuming that the vertical struc- 
ture of the gas follows a gaussian distribution, leads to a 
midplane density at 1 AU of 1.4 x 10^^ g cm^"^. 



The high density model: iDeschI (|2007l ) introduced a 're- 
vised MMSN model' by adopting the starting positions 
of the planets i n the 'Nice' model of planetary dynamics 
(iTsiea nis et al.] (j2005l )) thus taking into account planetary 
migration. The model predicts that the solar system started 
out in a much more compact configuration and its surface 
density profile is given by: 



Eg(r) 5.1 X 10^ 



cm2 VAU/ 



-2.2 



(4) 



This model is consistent with a decretion disk which is 
being photoevaporated by the central star. Although the 
model of Desch (200|) was defined for the outer solar 
system, we extrapolate the profile to 1 AU in order to 
cover a broad range of surface density values in our 
calculations. Assuming, as in the MMSN model, a gaussian 
vertical distribution, the density at 1 AU in the midplane 
is 2.7 X 10"^ g cm"3. 

For simplicity, we adopt a midplane temperature of 200 
K (isothermal sound speed of Cs = 8.5 x 10"* cm s^^) in all 
the three models. 



2.2. Relative velocities 

We consider three sources for relative velocities between 
dust aggregates. These are Brownian motion, radial drift 
and turbulence. In the following, we discuss these sources. 

The average relative velocity of two particles with mass 
TOi and 7712 in a region of a disk with temperature T due 
to Brownian motion is 



A7;b(771i, m2) 



' 8fcr(TOi -i- TO2) 

7rmi7)7,2 



(5) 



For micron sized particles, the relative velocity is of the or- 
der of 0.1 cm s~^, but for cm sized particles this value drops 
several orders of magnitude. Therefore, Brownian motion is 
only effective for collisions between small particles during 
the initial stages of growth. Coagulation due to Brownian 
motion results i n fiuffy agg r egates of fractal dimensio n 
around 2 and 3 (iBlum et al.i (1 19961): iKempf et all (|l999l ): 
iBlum et all (|2000l ): iKrause fc BlumI (|2004^ ). In practice 
there is no growth due to Brownian motion for aggregates 
larger than 100 micron. 

The second source for relative velocity is turbu- 
lence. Relative velocity of aggregates due to the ran- 
dom moti on of turbu l ent ed di es were calcu lated numer- 
icallv bv IVolket aLl (Il980l ). iMizuno et aH (fl988 ) and 
iMarkiewicz et al.l l[1991). We use t he cl osed form expres- 
sions presented bv lOrmel fc Cuzzil (|2007). We assume that 
turbu lence is parameterized by the Ishakura fc SunvaevI 
(|l97i) a parameter 



9' 



(6) 



where vt is the turbulent viscosity, Cg is the isothermal 
sound speed and Hg is the pressure scale height of the disk. 
The value of the a parameter reflects the strength of the 
turbulence in the disk. Typical values of a in this paper 
range between 10~^ and 10^'''. The turbulent relative ve- 
locity is a function of the stopping times of the two colliding 
particles. The stopping time (or friction time) is the time 
the particle needs to react to the changes in the motion of 
the surrounding gas. As long as the radius of the particle 
is smaller than the mean free path of the gas (a < fAmfp), 
the par ticle is in the Epstein regime, where the stopping 
time is (|Epsteinl (|1924l )): 



3777 



Ep 



(7) 



where 771 and A are the mass and the cross section of the 
particle, Pg and Uth are the gas density and the thermal 
velocity. At high gas densities, where the mean free path 
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10"'° 10"^ 10"^ 10"* 10"^ 10° 10^ 10"'° 10"^ 10"^ 10"" 10"^ 10° 10^ 

mass [g] projectile mass [g] 



Fig. 1. The combined relative velocities caused by Brownian motion, radial drift and turbulence for fluffy particles 
(5* — 20) in the three disk models for equal sized particl es (a) and for differe nt sized particles with a mass ratio of 
100 (b). The solid Hne indicates the low density model of iBrauer et al.l ()2008a[ ). Physical parameters of the disk: the 
distance from the central star is 1 AU, temperature is 200 K, the density of the gas is 2.4 x 10~^^ g cm~'^, and the 
turbulence parameter, a = 10~^. The dotted line represents the MMSN model. The density is 1.4 x 10^^ g cm~^, the 
other parameters are the same. The dashed line corresponds to the high density disk. The gas density is 2.7 x 10~* g 
cm~^. 



is low or in case of larger particles, the first Stokes regime 
applies and the stopping time is 



3rn 4 a 

ts = tgt = X . 

4:VthPgA 9 Amfp 



(8) 



In the first Stokes regime the stopping time is indepen- 
dent of the particle-gas relative velocity as well as the 
gas density. This regime can be used as long as the par- 
ticle Reynolds number is smalle r than unity. The parti cle 
Reynolds number calculated as (jWeidenschillind (|l977al ')): 



_ 2aAvpg 



(9) 



where Avpg is the relative velocity between the particle 
and the gas, and ij is the gas viscosity. For particles outside 
the Epstein regime, it can be assumed that the systematic 
velocity (radial drift) dominates over the random veloci- 
ties (turbulence); therefore, Avpg ~ I'd, where vd is the 
drift velocity of the particle, defined in the next paragraph. 
The particle Reynolds number never exceeds unity in our 
simulations. Therefore, we do not include further Stokes 
regimes. 

Radial drift also leads to relative velocities between ag- 
gregates. Radial drift {vd) has two sources: drift of individ- 
ual particles (vd) and drift due to accretion processes of the 
gas (vda), thus the total radial drift velocity is vd = v^+Vda- 
T he radial drift of i ndividu al dust aggregates with mass m 
is (|Weidenschillingl (jilizi)) 



2v 



Vd 



N 



St + l/St' 



(10) 



where St is the Stokes number of the aggregate [St = tgil, 
where f2 is the orbital frequen cy) an d vn is the maximum 
radial drift velocity (Whipple. ()1972D ). 

The second part of the radial velocity is due to the accre- 
tion of the gas. This part of the radial velocity is calculated 



o 
in 



10' 
10° 
10"^ 
10"" 
10"^ 
10"^ 



High density model 

MMSN model 

Low density model. 



10" 



particle radius [cm] 



Fig. 2. The Stokes number as a function of the particle 
radius in the three models. The parameters of the dust 
for all of the models are the following: monomer radius is 
flo = 0.75 /im, material density is po = 2 g cm~^, and 
= 1. 



as follows (jKornet et al.l (|200ir )): 



Vda 



'^gas 



i + st^' 



(11) 



gas 



wh ere f gas is the accretion velocity of the 
(Ta keiichi fc Lh] (|2002f l'l. 

The relative velocity due to radial drift is then simply 
the difference between the radial velocity of particle 1 and 
particle 2. However, as the Stokes number of the aggregates 
is always smaller than 10~^ (see Sec. S]), the second term 
of the radial velocity (vda) can be safely neglected: 



AvD 



\VD1 - VD2\ ~ \vdl - Vd2\- 



(12) 
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This study uses two quantities to describe the porosity 
of the aggregates. The volume filling factor is: 

<l^ = V*/Vtot^{A*/Af'\ (13) 

where V* is the volume occupied by the monomers and 
Vtot is the total volume of the aggregate, including pores, 
and A and A* are the surface area equivalents of these 
quantities. In this way, the filling factors also enters the 
definition of the friction time (Eqs. [7] and [8|). The density 
of aggregates then follows as p = po^j where = 2 g 
cm~^ is the material density of the silicate. In this study 
we will also use the reciprocal parameter of the filling fac- 
tor, which is denoted the enlargement parameter, 'J = 0^^. 

We illustrate the relative velocity between equal sized 
and different sized aggregates with \[' = 20 (0 = 0.05) in 
Fig.[T]for the disk models considered in this work. Adopting 
a threshold (fragmentation) velocity of 1 m s^^, the maxi- 
mum particle size, which can be reached in the models are: 
0.025 cm in the low density model, 1.4 cm in the MMSN 
model and 1.7 cm in the Desch model. The Stokes num- 
bers of these particles are the same in all the three models, 
4.7 X 10~^. The constant fragmentation velocity of 1 m s~^ 
is the typical velocity at which silicate particles will frag- 
ment (Birnstiel et al. (2009)). In our collision model this is 
not the case for all combinations of mass ratio and porosity 
(Paper I), but the m s~^ threshold is still a useful proxy 
for the point where fragmentation processes will become 
important. 

Figure [2] shows the Stokes number as a function of par- 
ticle radii in the three models. Initially, particles are in the 
Epstein regime, where the stopping time, thus the Stokes 
number, depends on the gas density. When the particles 
enter the Stokes regime, the stopping time becomes inde- 
pendent of the gas density (see Eq. [5]). One can see that 
particles in the Desch model are in the Stokes regime at 
Stokes number of 4.7 x 10"'^ (when the particles have rel- 
ative velocities of 1 m s~^), while the aggregates in the 
MMSN model are close to it, which explains why the max- 
imum particle size is almost the same in these two models. 

As discussed in lOrmel fc Cuzzil (|2007[ ). particles are ini- 
tially in the 'tightly coupled particle' regime, where the ed- 
dies are all of class I type, meaning that the turnover time 
of al l eddies is l onger than the friction time of the parti- 
cles (|V61k et al.l (jl980D ). A particle, upon entering a class I 
eddy, will therefore forget its initial motion and align itself 
to the gas motions of the eddy before the eddy decays or the 
particle leaves it. This regime is apparent in Fig. [T^ and b. 
Different sized particles are in this relative velocity regime 
as long as their masses are less than 10~® g in the low den- 
sity model, 10"^ g in the MMSN model and 10"^ g in the 
high density model assuming fluffy particles = 20). If 
the particles leave this regime and enter the 'intermediate 
particle' regime, their relative velocity increases. This tran- 
sition affects the particle evolution, as discussed in e.g. Sect 

3. Collision model and implementation 

In this work we use a statistical or 'particle in a box' method 
to compute the coUisional evolution. That is, we assume 



that all particles are homogeneously distributed within a 
certain volume (the simulation volume). In reality however, 
the particles could leave the simulated volume or new parti- 
cles could enter from outside due to radial drift or random 
motions (turbulence and Brownian motion). Since we do 
not resolve the spatial dependence of the aggregates, we 
will simply assume that local conditions hold during the 
run. The gas and dust densities are kept constant and par- 
ticles cannot leave or enter the simulation volume (hereafter 
'local approach'). 

3.1. Short overview of the collision model 

Many laboratory experiments on dust agg regate collisions 
have been performed in the past years, see iBlum fc WurmI 
( 2008 .) . The g rowth begins as fra ctal growth and we use 
the recipe of lOrmel et al.l (|2007l ) to describe this initial 
stage. However, once aggregates have restructured into 
non- fractal, macroscopic aggregates (e.g. > 100 /xm), 
laboratory experiments show that the coUisional outcomes 
become very diverse. In this regime, many new experiments 
were performed with dust aggregates consisting of 1.5 
fim diamete r Si02 monomers e ither with high porosity 
(f) = 0.15 ()Blum fc Schrapleil ()2004l )). or intermediate 
porosity (0 = 0.35). Paper I compiled 19 experiments 
with different aggregate masses, collision velocities, and 
aggregate porosities. 

From these experiments we have identified nine different 
coUisional outcomes involving sticking, bouncing, or frag- 
mentation (see Fig. The occurrence of these regimes 
mainly depends on aggregate masses and collision veloci- 
ties. However, it also depends on the porosity of the par- 
ticles and on the critical mass ratio. For example. Paper I 
finds fragmentation in collisions between a porous aggre- 
gate and a solid wall, whereas Langkowski et al. (2008) 
find sticking of a porous projectile by penetrating an also 
porous target. Likewise, Heifielmann et al., in prep, find 
bouncing of two similar-sized, porous dust aggregates, while 
Langkowski et al. (2008) find sticking for the same veloc- 
ity where one collision partner (target) was significantly 
bigger. To address the importance of the mass ratio and 
porosity, we have identified eight different collision regimes 
(look-up tables) based on a binary treatment of porosity 
and mass ratio: i.e., (i) similarly sized or differently sized 
collision partners and (ii) porous or compact collision part- 
ners. The further distinction between target, which we al- 
ways define as the heavier collision partner, and projectile 
then results in eight different collision regimes. We denote 
these regimes as 'pP' (porous projectile, porous target; tar- 
get significantly bigger than the projectile), 'pc' (porous 
projectile, compact target; target of similar size than the 
projectile), etc. 

In Paper I we have classified each of these 19 experi- 
ments into one or more of these eight regimes (see Fig. 10 
of Paper I). Based on extrapolation of experimental find- 
ings, we decide in which mass and velocity range collisions 
result in sticking, bouncing, or fragmentation. These results 
are presented in Fig. 11 in Paper I. 

It should be noted that the critical mass ratio be- 
tween the equal-size {'pp', 'cc \ etc.) and the different size 
regimes { 'pP', 'cC\ etc.) is ill-constrained by experiments. 
Therefore, we use critical mass ratios of r,„ = 10, 100 and 
1000 to explore the effect of this parameter. 



6 



A. Zsom et al.: The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals? II. 




B1 (bouncing with compaction) 



■O O- 



■o 




51 (hit & sticl<) ^ — V. 

u 

52 (sticl<ing through surface effects) 

o 



00 



S3 (sticlfing by penetration) 



S4 (mass transfer) 




o „= 




B2 (bouncing with mass transfer) 



F1 (fragmentation) 




o ^ 
o 



- 'o 
O n .O 



F2 (erosion) 



F3 (fragmentation with mass transfer) 




Fig. 3. The collision types considered in this paper. We distinguish between similar sized and different sized particles. 
Some of the collision types only happens for one of the mass ratios. Grey color indicates that during the given collision 
type the particle is compact, or part of the mass will be compacted. 



3.2. Porosity 

Paper I defined a binary representation of the porosity, par- 
ticles are either porous or com pact. Following the simple 
model of 'Weidh ng eTal] (1200^ 1. we include a continuous 
transition between these two 'phases'. They showed that 
the compaction of particles due to bouncing can be de- 
scribed by porous and compacted sites on the surface of 
the aggregate. A site of the aggregate is porous if it did not 
encounter any collisions yet (e.g. bouncing), a compacted 
site encountered at least one coUision already but any fur- 
ther coUision happening at that part of the surface cannot 
change the porosity of this site anymore. We describe the 
probability of hitting a passive site of the aggregate in the 
following way: 

P. = 4^, (14) 



leads to more compact stru ctures with f r actal dimension 
of 3. The hit&stick recipe of lOrmel et J] (j2n07D attempts 
to interpolate between these two fractal models. At one 
point, coUisional energies become large enough to invali- 
date this assumption. This occurs when the collision energy 
i s five times higher t han t h e rolling energy of mo nomers 
(iDominik fc TielensI (|1997r i. iBlum fc Wurn] (120001 ')'). The 
internal structure then becomes homogenous. In our model 
we assume that once the fractal growth due to SI is over, 
Bouncing with compaction (Bl) restructures the aggregates 
producing compact structures with fractal dimension of 3. 
As our model can not follow the exact shape of the par- 
ticles, we assume that the aggregates are spheres and can 
be described with a single density (p ) thus negl e cting the 
effects of e.g. the 'toothing radius' of lOssenkopj ()1993l ) or 
the craters forming during Penetration (S3). 



Using the expressions for the relative velocity, the col- 
lisional cross section between the dust particles, and 
the collisional outcome, we solve for the temporal evo- 
lution of the dust size distribution. Traditionally, the 
Smoluchowski equation is solved to follow the evolution of 



where (j) is the volume filling factor of the aggregate, is 3,3, The Monte Carlo method 
the critical porosity (0c = 0.4, see Paper I), and 0p is the 
volume filling factor of the porous site, which is chosen to 
be 0.15. If (j) is between 0.15 and 0.4, a random number 
decides whether the particle collided with a porous or a 
compact site. Such a treatment of the porosity ensures 
a continuous transition from porous to compact aggregates. 

During the initial Hit & Stick (SI) phase, particles are 
in the f ractal growth regime (lOssenkopj ()l993l ) , iBlum et"ari 
(i2000l ). lKrause fc BlumI (|2004l )). Particles grow initially due 
to Brownian motion and later due to turbulence. The struc- 
ture of the aggregate depends on whether the collision 
happened between a cluster and a monomer (PCA) or 
between two clusters (CCA). The latter results in fiuffy 
aggregates with fractal dimension of 2, while the former 



the m ass distribution function (e.g. DuUemond & Dominik 
(|2004r) . DuUemond & Dominik (2005), Tanaka et al 



20051 ). iBrauer et al.l (|2008aD ). The continuous form of 



the Smoluchowski equation used in these wo rks lacks 
the s tochasticity of the coagulation problem (|Safronovl 
(|l969f )). All bodies with mass m will grow in the same way 
thus the spatial and temporal fluctuations of the particle 
ensemble are averaged out. In reality, however, particles 
with similar masses might follow a different evolutionary 
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path depending on what other particles they collide 
with. The collision model typically used in these works 
is, by necessity, rather simple as in the Smoluchowski 
formulation the collision and time evolution steps are 
linked together. These collision models consist of sticking 
and fragmentation and only the mass of the particles 
is followed. The advantage of such a model is that it is 
computationally not too expensive: The entire disk can be 

mo deled. 

lOrmel et al.l ()2007l l introduced a new Monte Carlo 
method to solve for the mass and the porosity distribu- 
tion function simultaneously. Their collision model consists 
of sticking and compaction; ZsD08 added a simple frag- 
mentation model as well. Although these models are more 
detailed, one can see that they still lack the full complex- 
ity which is observed at "the zoo" of laboratory collision 
experiments. 

The MC-approach used in this study has previously 
been presented by ZsD08. It can be characterized by two 
key properties: (1) the number of MC-particles (also re- 
ferred to as representative particles) is kept constant; (2) 
the method follows the mass of the particle distribution. 

Property (1) is required to preserve good statistics. 
Because of the ^fN noise of MC-methods, a large fluctua- 
tion of N would s e verely affect the accuracy of the method 
(|Ormel fc Spaansl ()2008l )). The second property states that 
our primary interest lies in the particles that contain most 
of the mass of the system. Moreover, it has been shown that 
following the particle's mass distribution - rather than the 
number distribution - is also a prerequisite to preserve a 
good co rrespondence with syste ms that experience strong 
growth (jOrmel fc Spaansl (|2008l )). 

Property (2) ensures that the MC method samples the 
parameter space only where a signicant portion of the total 
dust mass is. However, this is not always desirable. For 
instance, radiative transfer calculations require the surface 
area distribution of the aggregates, which determines the 
opacity. If most of the particle mass is contained in big 
particles (which are not observable) the amount of small 
particles (which could contain most of the surface area and 
determines the IR appearance of the disk) might be resolved 
with a bad statistics. But if we are interested in following 
the evolution of the dominant portion of the dust, then MC 
methods naturally focus on these parts of the phase space. 

A required condition for the ZsD08 method to work is 
that the number of the representative particles N is much 
less than the number of actual aggregates present in the 
system under consideration - a condition that is safely 
met in any of our simulation runs. Then, a representative 
particle will collide only with the non-representative 
particles, whose distribution is assumed to be the same as 
that of the representative particles. We refer to ZsDOS for 
details regarding the precise implementation and accuracy 
of the method; here we further concentrate on how the 
method operates under the new coUisional setup. 

The collision kernel is defined as the product of the cross 
section of the colliding particles and their relative velocity: 

^i,fc = cTi.fcAwj^fc, (15) 

where the index i corresponds to the representative particle 
and k is the index of the non-representative particle. The 
kernel is proportional to the probability of a collision. The 



value of i^i.fc is calculated for every possible particle pair, 
and random numbers determine which of the collision will 
occur first and at which time interval. 

The above properties and conditions specify the essence 
of the ZsDOS method: one of the two collision particles is a 
representative particle and, by property (1), only one of the 
coUisional products becomes the new representative parti- 
cle. By property (2) the choice for the new representative 
particle is weighed by the mass of the collision products. 
A very helpful analogy here is that of the representative 
'atom', which is contained within the representative par- 
ticle. The choice for the new representative particle after 
the collision is then proportional to the probability of the 
representative 'atom' ending up in the collision products. 
If, for instance, a collision leads to the production of an en- 
tire distribution of debris particles, the probability that a 
particular debris fragment becomes the new representative 
particle is proportional to the likelihood of this fragment to 
contain the representative 'atom'. 

3.4. Implementation of the collision types 

We describe the implementation of the collision model using 
the representative 'atom' concept. We refer to Paper I for 
details of the various collision types described below. 

Hit & Stick (SI), Sticking through surface effects (S2), 
Penetration (S3): All three of these collision types result 
in sticking and increase the mass of the aggregate by that 
of the projectile, but the porosity changes in a different 
manner (see Paper I) . The new mass of the representative 
particle i is then the sum of the original particle masses, 
Winnow = -|- mfe, where rrii is the mass of the representa- 
tive particle and is the mass of the non-representative 
particle. 

IVlass transfer (S4): In the case of Mass transfer (S4), a 
certain percentage of the mass of the projectile sticks to 
the target, while the left-over mass of the projectile will 
fragment into a power law distribution (see Paper I). 
There are two situations to consider: 

1. The representative 'atom' is part of the target. The 
mass of the new aggregate will be the mass of the orig- 
inal aggregate plus the transferred mass from the non- 
representative particle (m^^new = + "i-ti-ans: whcrc 
'Titrans is the transferred mass calculated according to 
Paper I). 

2. The representative 'atom' is part of the projectile. 
Again, we have two situations. 

(a) The representative 'atom' will be transferred to the 
non-representative particle. The mass of the new 
representative particle will be the mass of the non- 
representative particle plus the transferred material 
('7ii,now = ™/c + Wtrans)- The probability of trans- 
ferring (removing) the representative atom from the 
projectile is simply P = mtrans/wz, the ratio be- 
tween the transferred mass and the mass of the pro- 
jectile. 

(b) The representative 'atom' remains in one of the 
fragments. The probability of this event is P = 
[m,i — mtrans)/wi, the ratio between the fragmented 
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mass to the original mass of the representative par- 
ticle. As discussed in Paper I, the fragments follow a 
power law mass distribution. The distribution is de- 
fined by the maximum mass of the fragments, which 
is a function of the relative velocity and the total 
mass of the fragments. The total mass of the frag- 
ments is rrii — mtrans- We randomly choose from the 
fragment distribution to find the new mass of the 
representative particle (to find which of the frag- 
ments will contain the representative 'atom'). 

Bouncing with compaction (Bl): Upon Bouncing with com- 
paction (Bl) particles collide and bounce. Bouncing itself 
does not change the mass of the particles, but it compact- 
ifies them ac cording to Pape r I. As observed in laboratory 
experiments (jWeidling et al.l (|2009l )). there is a small prob- 
ability (Pfrag = 10~^) that the bouncing particle will break 
apart. If this happens, we break the particle into two equal 
mass pieces. 

Bouncing with mass transfer (B2): Bouncing with mass 
transfer (B2) is, from the implementation point of view, 
similar to Mass transfer (S4). The recipe to define the 
new representative particle is as in Mass transfer (S4). The 
difference is that the projectile does not fragment during 
the collision, and that the porosity changes differently (see 
Paper I). 

Fragmentation (Fl): Fragmentation only happens between 
similar sized aggregates in the 'pp' and 'cc' regimes. The 
fragments follow a power law mass distribution where the 
maximum mass of the fragments is determined by the rel- 
ative velocity of the particles and the total mass that goes 
into the fragments (Paper I). We randomly choose from 
these distribution to determine the new mass of the repre- 
sentative particle. 

Erosion (F2): Erosion (F2) happens between different sized 
particles only. During the collision the projectile "kicks out" 
pieces from the target aggregate. These pieces follow a 
power law distribution (see Paper I). We have to consider 
two cases. 

1. The representative 'atom' is in the target. Again, we 
have two possibilities. 

(a) The representative 'atom' will stay in the target af- 
ter the collision. The mass of the new particle will 
be rrii^new = mi — mer, where rucr is the eroded mass. 
The probability of this event is P = (mj — mf.,:)/mi, 
that is the ratio between the left-over mass (which 
does not erode) and the mass of the original particle. 

(b) The representative 'atom' is part of the eroded par- 
ticles. As the eroded particles follow a power law dis- 
tribution, we randomly pick from this distribution to 
determine the new mass of the representative parti- 
cle. The likelihood for this event is the ratio between 
the eroded mass and the original mass of the particle 
(P = mcr/mi). 

2. The representative 'atom' is part of the small particle 
which caused the erosion. As the particles do not stick 
and the small particle does not fragment, the represen- 
tative particle remains unaffected. 



Fragmentation with mass transfer (F3): In Fragmentation 
with mass transfer (F3) the porous particle gets destroyed 
by the compact one and transfers a certain amount of mass 
to the compact particle. Fragmentation with mass transfer 
(F3) only happens in the 'cp ' regime. Again, we have two 
possibilities. 

1. The representative 'atom' is part of the compact parti- 
cle. In this case the representative 'atom' cannot leave 
the particle. The new mass of the representative parti- 
cle will be mi. new = m,j + mtrans: the sum of the original 
mass plus the transferred mass. 

2. The representative 'atom' was part of the porous aggre- 
gate. 

(a) The representative 'atom' is part of the material 
which is transferred to the compact particle. In this 
case, the new mass of the particle will be that of 
the compact (non-representative) particle plus the 
transferred material (m^^ncw — "m-k + mtrans)- The 
probability of this event is P = mtrans /m^. 

(b) The representative 'atom' is part of the fragments. 
As before, the mass distribution will follow a power 
law and we randomly pick from this distribution to 
determine the new mass of the representative par- 
ticle. The probability of this event is P = (m^ — 

mtrans) /mi. 

3.5. Evolving the particle properties in time 

We summarize how the particle properties are evolved in 
time using the above described kernel. We start with the 
size and porosity distribution of the particles at a given 
time, t. At t = 0, we must give the initial size and porosity 
distribution, see Sect. 14.1] Knowing these: 

— We calculate the cross sections of all possible collision 
partners, as well as their relative velocities using the 
equations described in Sect. 12.21 Both are used to de- 
termine the collision rates between the particle pairs. 

— By using random numbers, we identify from the collision 
rates the representative particle, which is involved in 
the collision, as well as the non-representative particle 
it collides with and at what time the collision takes place 
(t + At). 

— Knowing the masses (mass ratio) and porosities of the 
collision partners, we identify in which of the eight 
regimes the collision takes place (e.g. 'pP', or 'pC, 
etc.). 

— Next, we identify which of the nine collision types ma- 
terializes (Fig. |3| using the relative velocity of the par- 
ticles and the mass of the projectile (see Paper I). 

— Based on the collision recipe described in Paper I and 
Sect. 13.41 the new mass and new porosity of the rep- 
resentative particle is calculated and the new size and 
porosity distribution of the particles at time t 4- At is 
obtained. 

— In the final step, we update the collision rates. 

3.6. Numerical issues 

As mentioned in ZsD08, a sufficiently high number of rep- 
resentative particles is needed to properly reproduce the 
physics of the collision kernel. We performed simulations 
with an increasing number of representative particles and 
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found that for more than 200-300 particles the results 
of the simulations do not change significantly. For all of 
the simulations described in the following sections 500 
representative particles are used and we average the results 
of 20 simulations to decrease the numerical noise of the 
code. 

The required computational time strongly depends on 
the collision rate of the particles thus determined by the 
dust density and the relative velocity (the a turbulence 
parameter mostly) and the length of the simulation. On a 
2.83 GHz CPU, running the simulations on a single core, the 
CPU time varies between nine hours (the low density model 
with a = 10^^) and three days (the high density model with 
a — 10^"^). Both simulations covered 10^ years of particle 
evolution. The high density simulation with a = 10"'*, crit- 
ical mass ratio of 100 and t = 10^ years of evolution takes 
twelve days to simulate. 

4. Results 

4.1. Initial conditions, setup of simulations 

All simulations start with silicate monomers of 1.5 fim di- 
ameter and 2 g cm~^ material density (monodisperse size 
distribution). We simulate the dust evolution at the mid- 
plane of our disk models at a distance of 1 AU from the 
central star. The gas density is obtained from the disk mod- 
els described in Sect. 12.11 We assume a typical 1:100 dust 
to gas ratio. We follow the history of each collision: the 
mass and porosity of the colliding particles, their relative 
velocity, the occurred collision type and the new mass and 
porosity of the particles. In this way we can reconstruct the 
history of the dust evolution. 

The parameters we vary in this study are the gas density 
Pg and the turbulence parameter a. We also treat the crit- 
ical mass ratio as a free parameter in order to explore 
its effect on the dust evolution. 

We provide a detailed description of the low density 
model with a = 10""* and critical mass ratio of 100 in Sect. 
14.21 We then compare this with the MMSN model and the 
high density model using the same turbulence parameter 
and the critical mass ratio (Sects. 14.31 and 14. 4p . In Sects. 
14.51 and 14.61 we discuss the effects of changing the turbu- 
lence parameter and critical mass ratio by comparing those 
results with the two example runs. 

4.2. The low density model 

The gas density in this disk model at 1 AU is 2.4 x 10~*^ 
g cm~^, the turbulence parameter is a = 10^*, the critical 
mass ratio is = 100. As shown in Fig. [U the particles 
reach the fragmentation velocity (1 m s~*) already at sizes 
smaller than millimeter because the particles in low gas 
density environment decouple from the gas already at these 
small radii. 

Figure |3] shows the evolution of the mass distribution 
(a) , the porosity distribution (b) and the collision frequency 
of the various collision types (c) . The x-axis shows the time 
in a logarithmic scale. The y-axis of Fig. 2^, b shows the 
mass and enlargement distributions, respectively. Here, the 
intensity of the color reflects the number density of repre- 
sentative particles, which, as explained in ZsD08, measures 
the mass density of the distribution. Thus, in Fig. the 



intensity levels directly reflect the mass density, while in 
Fig. EJd the colours indicate the mass weighted enlargement 
parameter. The black lines show the average of these quan- 
tities over the particle distribution. The y-axis in Fig. HI; 
represents the nine collision types used in this paper. Every 
stripe shows the total collision rate of the collision types at 
a given time. 

Figure [5] represents the collision history in the eight col- 
lision regimes. The x-axis is the velocity, the y-axis shows 
the mass of the projectile. A mass- velocity grid is created 
and for all grid cells we calculate how many collisions hap- 
pened inside that given grid cell during which our 'local 
approach' assumption is correct, that is 4 x 10^ yr. The 
different collision types and their border lines, as well as 
the areas which are covered with laboratory experiments 
(indicated with grey colors) are plotted. For more details 
on the experiments, see Paper I. 

In the 'cc' panel, we indicate two curves with dotted 
lines. These curves are evolution tracks. The left curve is 
obtained by calculating the relative velocity between equal 
sized particles with an enlargement parameter of 2.5 (vol- 
ume flUing factor of 40%). The right curve represents the 
relative velocity between particles having a mass ratio of 
100. These two curves serve as a guide to our results, as 
collisions should happen between these two curves in the 
'cc ' panel. The lower part of the left curve, where the rela- 
tive velocity decreases with increasing mass, is a sign that 
relative velocities between equal sized particles are domi- 
nated by Brownian motion. For higher masses, the relative 
velocity is dominated by turbulence. These curves do not 
precisely match the contours because we assumed a con- 
stant enlargement parameter of 40% when calculating the 
evolution tracks, whereas ^' is a free parameter in the sim- 
ulation. 

4.2.1. Early evolution 

We discuss here the evolution of the distribution functions 
until the 'local approach' assumption becomes invalid (4 x 
10^ yr). The long term evolution of the dust is discussed in 
Section 

We distinguish two distinct phases here. During the 
first 300 yr, particles grow by the Hit & Stick (SI) mecha- 
nism. The second phase is Bouncing with compaction (Bl) 
dominated; the particles leave the SI regimes. During this 
phase the mass of the particles is slowly decreasing and the 
enlargement parameter asymptotically reaches a minimum 
value of 2.23. As discussed in Paper I, keeping the bouncing 
velocity of a particle constant, the porosity of the aggre- 
gate will asymptotically reach a maximum value, 0max (see 
Paper I) . The relative velocity of a particle is a function of 
the friction time (Eq. [T]) , which depends on the ratio of the 
mass to surface area, m/A. Since particle growth is halted 
at this point in the simulation (m stays constant), only a 
decrease in A due to compaction can further increase the 
velocity between particles. The particle radius can decrease 
until either (/)max for the given relative velocity is reached, 
or until particles reach the maximum compaction possible. 
The latter limit, random close packing (RCP), corresponds 
to an enlargement parameter of 1.6 (volume filling factor 
of -60%). 

We find that fragmentation does not play a role during 
the evolution of these particles indicated by Fig. As 
can be seen in Fig. [51 their evolution is halted by bouncing 
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Fr. with mass transfer - F3 H 
Erosion - 1^2 H 
Fragmentation - F1 H 
B. with moss tronsfer - B2 H 
Bouncing with compaction - B1 H 
Mass transfer - S4 H 
Penetration - S3 H 
St. through surface effects - S2 H 
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Fig. 4. The evolution of the mass distribution (a), enlargement parameter distribution (b) and the collision frequency of 
the nine different collision types (c) in the low density model with a = 10~^ and critical mass ratio of 100. The x-axis 
shows the time. The y-axis of the (a) and (b) figures show the logarithmic mass and the linear enlargement parameter 
respectively. The contours represent the normalized mass density and the mass weighted enlargement parameter. The 
black lines represents the average of the mass and enlargement parameter at a given time. The y-axis on the (c) figure 
represents the nine collision types. Each stripe shows the total collision rate of the collision types. Two distinct phases 
can be distinguished. During the initial 300 yr particles grow by Hit & Stick (SI), after that the evolution is governed by 
Bouncing with compaction (Bl). The white lines indicate how long our 'local approach' assumption is valid (discussed 
in Sect. 131). 
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Fig. 5. The collision history of the eight regimes in the low density model for a — 10~^. The x-axis is the relative velocity, 
the y-axis shows the projectile mass. The different collision types, their border lines, as well as the areas covered with 
laboratory experiments (grey) are plotted. A relative velocity - mass grid is created and in these grid cells we calculate 
how many collisions happened until the 'local approach' assumption is valid (4 x 10^ yr). This is represented by the colors: 
yellow and red indicate a high collision frequency. The two dotted lines on the 'cc' regime are evolution tracks. Assuming 
a constant (40%) volume filling factor, the relative velocity between equal sized particles (left curve) and particles with 
a mass ratio of 100 (right curve) can be calculated. The collisions in the simulation should lay between these two lines. 
The small deviations are due to the fact that the volume filling factor is not exactly 40% during the simulation. 
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before the particles could reach the fragmentation barrier. 
The two dominant collision types are Hit & Stick (SI) and 
Bouncing with compaction (Bl). 



4.2.2. Termination of growth 

As we can see from Fig. [3 sticking at higher energies than 
the Hit & Stick (SI) border lines is only possible inside the 
'pP' regime. As soon as we no longer have collisions inside 
this regime or the SI regimes, the growth is halted. There 
can be two reasons why this is happening: 1.) All parti- 
cles are compact; there are simply no collisions in the 'pP' 
regime. 2.) The width of the particle mass distribution is 
less than the critical mass ratio (?"„), such that all collision 
take place in the equal-size regimes {'pp', 'pc\ etc.). 

In the case of the current simulation, the small particles 
have been 'consumed'. Once the heavy particles grow into 
the Bouncing with compaction (Bl) area of the 'pp ' regime, 
their growth in the 'pp' regime stops. The heavy particles 
collect the small ones via collisions in the 'pP' regime and 
by doing so, the width of the distribution is reduced to 
a value which is less than r^. Therefore, before particles 
could reach the fragmentation barrier, growth is halted. 
Due to Bl, particles get compacted and collisions in the 
'cc', 'cp' and 'pc' regimes appear. 

4.2.3. Long term evolution 

Before discussing the long term evolution of the distribution 
functions, we must consider for how long our starting as- 
sumptions ('local approach' and constant gas density) hold 
true. 

Using Eq. [TlJ we calculate that a particle with Stokes 
number 10^* drifts a distance of 1 AU in roughly 4 x 10^ 
yr. This is the drift timescale beyond which the 'local ap- 
proach' assumption (discussed in Sect.[3|) is not valid any- 
more: particles become separated from each other on this 
timescale. 

Another process through which particles separate is by 
viscous spreading. We determine the viscous timescale of 
the disk at 1 AU: 

Uis = r'lvT. (16) 

where r is the distance from the central star (1 AU), vt is 
defined in Eq. El The viscous timescale in our model, using 
a = 10"'*, is of the order of 10^ yr. 

One has to consider the results of the simulation 
with caution for longer times than the drift or viscous 
timescales. The equilibrium or final state of the particles 
is reached when mass decrease during Bouncing with 
compaction (Bl) and Bouncing with mass transfer (B2) 
and mass increase by Hit & Stick (SI) and Sticking through 
surface effects (S2) are in equilibrium. Or in other words, 
the final state is reached when the evolution of the average 
mass and the enlargement parameter is only determined 
by the stochastic fluctuations of the simulation. We find 
that the equilibrium state of the particles is hardly reached 
within these timescales. Upon neglecting these warnings, 
we find that the equilibrium state of the dust is reached 
at t = 4 X 10^ yr. The equilibrium is reached between 
the bouncing collisions resulting in breakage and Hit & 
Stick (SI) (see Fig. St). The equilibrium average mass and 



porosity of the particles are mgn = 2 x 10 * g, ^Pfin = 2.77. 

To be able to compare the distribution functions of dif- 
ferent runs, we define some quantities using the mean of 
the distribution functions shown with black lines in Fig. 2^ 
and b: max(m), the maximum of the mean mass; max(5'), 
maximum of the mean enlargement parameter; "^mim the 
minimum mean enlargement parameter when particles do 
not compact anymore; inoc, the time when ^'mm is reached, 
that is when the time derivative of ^ is zero [d^ jdt = 0); 
and max(5t), the maximum average Stokes number reached 
during the simulation. The values of these quantities are 
listed in Table [1] (model id 'Ltld-4ml00'). In this table, 
Col. 1 describes the model names. 'L' stands for the low 
density model, 'M' is the MMSN model, 'H' is the high den- 
sity model, the letter 't' and the following number indicates 
the value of the turbulence parameter, and the letter 'm' 
and the number shows the used critical mass ratio values. 
Columns 2, 3 and 4 show the gas density, turbulence pa- 
rameter and the critical mass ratio respectively. Columns 
5, 6, 7, 8 and 9 list the parameters defined to character- 
ize the distribution functions. These are max(m) in Col. 5, 
max(^') in Col. 6, '^min in Col. 7, t„oc in Col. 8 and finally 
max(5't) in Col. 9. 

4.3. The MMSN model 

The gas density in the MMSN model at 1 AU at the mid- 
plane is 1.4 X 10~^ g cm~^, a = 10~^, the critical mass ratio 
is 100. As shown in Fig.[Tl the particles grow to bigger sizes 
than in the low density model, as they are better coupled 
to the gas and the relative velocities are suppressed. As in 
the previous Section, we first discuss the evolution of the 
distribution functions for as long as the 'local approach' 
assumption holds true (6 x 10^ yr in this model). 

Figure [S] shows again the time evolution of the mass (a), 
enlargement parameter (b), and the collision frequency (c). 
Figure [7] shows the collision history. These figures show a 
rather different evolution than the previous model. 

4.3.1. Early evolution 

We find that during the fractal growth regime, the collision 
rate of Hit & Stick (SI) is much higher than in the low 
density model (Fig. ISt). This is due to the higher dust 
densities. We can see from Fig. [7l 'cc' regime, that growth 
starts with Brownian motion because the relative velocity 
decreases with increasing particle mass for particle masses 
less than 10~^ g. As a result of these low velocity collisions, 
some particles reach enlargement parameter values higher 
than 30 (volume filling factor less than 3.3%). At 200 yr, 
some particles grow above the border line of Hit & Stick 

(51) and enter the area of Sticking through surface effects 

(52) in the 'pP' plot, and Bouncing with compaction (Bl) 
in the 'pp' plot. Growth due to SI and S2 continues until 
different sized particles enter the transition regime in the 
'pP' plot. One can see in Fig. [B^, that some particles reach 
1 g in mass. However when particle collisions enter the 
transition regime between Bouncing with mass transfer 
(B2) and Sticking through surface effects (S2) in the 'pP' 
plot, their masses are equalized due to the mass transfer of 
the B2 collisions and the collisions shift to the similar sized 
regime (Bl). We find that after roughly 10^ yr particles 
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10° 10' 10^ 10^ 10" 10^ 10® 10' 
time [years] 



Fig. 6. Same as Fig. |5]but for the MMSN model. We magnify the spike of the mass distribution at ~ 2.5 x 10^ yr in 
Fig. (a). Four phases can be distinguished here. Initially (first 300 yr) particles grow purely by Hit & Stick (SI). After 
this the growth slows down because Bouncing with compation (Bl) starts and all particles leave the Hit & Stick (SI) 
regime. Between 3 x 10'' and 10^ yr, particles enter the transition regime between Sticking through surface effects (S2) 
and Bouncing with mass transfer (B2) on the 'pP' regime. Some particles reach masses of 1 g, but their masses are lastly 
reduced by B2. The last phase is Bouncing with compaction (Bl) dominated. The solid/dotted white lines indicate how 
long our 'local approach' assumptions are valid (discussed in Sect. [3]). 
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Fig. 7. Same as Fig. [5] but for the MMSN model. The particles are better coupled to the gas due to the higher 
density. Therefore, they grow to bigger sizes than in the low density model. 
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mostly bounce and compact. The enlargement parameter 
reaches a minimum value of 1.85 (54% volume filling 
factor), the mass distribution function slowly decreases 
due to a small probability of breakage. Collisions at this 
point are mainly happening in the 'cc' regime. 

A peculiar feature of Fig. [SK is a peak at t = 2.5 x 10"^ 
yr, which is accompanied by a fast decrease in the enlarge- 
ment parameter in Fig. [B)d and an increased collision rate 
of Sticking through surface effects (S2) and Bouncing with 
mass transfer (B2) in Fig. At this point, the relative 
velocity due to turbulence increases. As discussed in Sect. 
12.21 particles leave the 'tightly coupled particle' regime and 
enter the 'intermediate particle' regime (see the relative ve- 
locity bump in Fig. [TJ)). We calculate the growth timescale 
of the heaviest particle with mass M in the simulation as 
follows: 

/ 1 dM\~^ 

This is illustrated with a dotted line in Fig. [H As a com- 
parison, we also calculate the minimum growth timescale 
that a particle can have (solid line). That is: 



1000 F 



\pdAvaM J 



(18) 



where um is the cross section of the largest particle. Here, 
we assume that the 'swept up' particles have masses of 
M/lOO, therefore we use the relative velocity curve pre- 
sented in Fig. [TJ), dotted line. The effect of the relative 
velocity bump and the increased growth rate is seen at 0.1 
g- 

The relative velocity 'boost' happens shortly after the 
particles enter the transition regime of S2 and B2 in the 
'pP' plot. The heaviest particle, which encounters the 
velocity transition the earliest, experiences higher relative 
velocities leading to an increased collision rate with the 
other particles. As the particles are initially located at the 
lower part of the S2-B2 transition regime (with masses 
of 10~^ g, see Fig. |B^), the heaviest particle experiences 
fast growth and reaches masses of 30 g. The simulated 
timescale, however, does not reach the minimum growth 
timescale due to the Bouncing with mass transfer (B2) 
collisions which are reducing the mass of the heaviest 
particle. The rest of the particle population increases in 
mass because of B2 and the growth rate of the heaviest 
particle decreases. Eventually, the fast growth of the 
heaviest particle is halted, the growth timescale at m = 30 
g is infinity. From this point on, the heaviest particle re- 
duced in mass, and B2 equalizes the masses of the particles. 



4.3.2. Long term evolution 

We calculate the drift and viscous timescales to determine 
how long our assumptions of 'local approach' and constant 
gas density are valid. Assuming Stokes number 10~* par- 
ticles, we find that the drift timescale is of the order of 
6 X 10^ yr, the viscous timescale is 10^ yr. These timescales 
are indicated with solid and dotted white lines in Fig. [SI 

We find that the final equilibrium is reached at 
t — 2 X 10^ yr, which is longer than the drift and the 
viscous timescales. The equilibrium is reached between the 



100 



o 

"D 




Fig. 8. The dotted line and the '-I-' signs represent the 
growth timescale of the heaviest particle in the MMSN sim- 
ulation with a = 10^"* and — 100. As a comparison, we 
show the minimum growth timescale a particle can have in 
this simulation (solid line). 



growth mechanisms of Hit & Stick (SI), Sticking through 
surface effects (S2) and the destruction mechanisms of 
bouncing resulting in breakage and Bouncing with mass 
transfer (B2). The final average jnass and porosity of the 
particles are man = 2 x lO"-^ g, ^fin = 3.3. 

We conclude that the dust evolution is more complex in 
the MMSN model than in the low density model because 
the complex interaction of the velocity field and the col- 
lision kernel is apparent in this model. As in the previous 
model. Bouncing with compaction (Bl) is the most frequent 
collision type and Hit & Stick (SI) determines the initial 
particle growth, but Sticking through surface effects (S2) 
and Bouncing with mass transfer (B2) are of importance in 
this model. The final equilibrium is not reached within the 
drift and viscous timescales. 



4.4. The high density model 

The gas density in this model is 2.7 x 10~^ g cm~^ at the 
midplane of the disk at 1 AU distance from the central star. 
The values of a, and the dust to gas ratio are the same 
as in the previous models. 

Figure [U dashed line, shows the relative velocity field 
of fluffy aggregates in this model. As already discussed in 
Sect. 12.21 the aggregates reach 1 m s~^ relative velocities 
at similar masses as the MMSN model due to the Stokes 
drag. Therefore, we expect that the final aggregate sizes 
and masses will be similar to the particles produced in the 
MMSN model. 

Figure |9] shows the time evolution of the mass (a) , en- 
largement parameter (b) and the collision frequency (c). 
Figure [TUl illustrates the collision history. 

4.4.1. Early evolution 

As seen in Fig.llO( Brownian motion is the dominant source 
of relative velocity, as long as particles stay below masses 
of 10~* g (that is an order of magnitude higher than in 
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the MMSN model). Therefore, the enlargement parameter 
of the aggregates is also higher than in the MMSN model. 
As the Hit & Stick (SI) collisions are more frequent than 
in the MMSN model due to the higher dust densities, the 
particles reach the Sticking through surface effects (S2) - 
Bouncing with mass transfer (B2) transition regime earlier, 
at t = 200 yr. The peak in the mass distribution is not as 
pronounced as in the MMSN model. The relative velocity 
boost happens for heavier aggregates (10~^ g, see Fig. [T}d) 
due to the higher gas density of the model. When the fast 
growth of the heaviest particle starts, most of the projectiles 
are already in the transition regime. Here, the B2 collisions 
soon reduce the mass of the heaviest particle and narrow 
the mass distribution. 

In contrast to the MMSN model, the mass of the parti- 
cles is not reduced due to the low probability of breakage in 
Bouncing with compaction (Bl), but is kept nearly constant 
in time. This is the result of the increased collision rate of 
Sticking through surface effects (S2). The S2 collision rate 
increased because of low velocity collisions, which are oc- 
curring when particles are in the tightly couple regime and 
have similar stopping times. These S2 collisions are happen- 
ing in the 'pp' regime as seen in Fig. [TU] These collisions 
cancel out the effect of breakage in Bl. 

The maximum Stokes number reached in this model is 
3.6x 10~5 (see TableH model id 'Htld-4ml00'), lower than 
in the MMSN model. The growth in this model is halted by 
the Bouncing with mass transfer (B2) collisions in the tran- 
sition regime of the 'pP' panel. This shows us that particles 
cannot reach masses much larger than 1 g independently 
from the gas density (or Stokes number), because at this 
point, particles enter the S2-B2 transition regime and the 
growth is halted. Further increasing the gas density would 
result in even lower Stokes numbers. 

4.4.2. Long term evolution 

The drift and the viscous timescales in the high density 
model are both 10^ yr. As seen in Fig. [3^, the particle 
masses do not change significantly after t = lO'^ yr. The 
porosity is reduced due to Bouncing with compaction (Bl) 
and it reaches a final value of 5.41 a,t t — 3 x 10^ yr. 

4.5. Varying the turbulence parameter 

To explore the effects of turbulence, we perform two more 
simulations in each of the disk models. We keep the critical 
mass ratio fixed (100) and vary only the turbulence param- 
eter (a) to have values of 10"'^, lO"'' and 10^^. The results 
are shown in Table [1] the first nine models and in Fig. [TTk . 

The work of lBrauer et al.l (|2008al ) suggests that in situa- 
tions where fragmentation limits the growth, a lower turbu- 
lence strength results in bigger aggregates. This, of course, 
directly reflects the shift of the fragmentation threshold (1 
m/s) to larges sizes when a is lower (Fig. [T]). In this study 
it is fragmentation that balances the growth, which results 
in a (quasi) steady-state. For the low density models we do 
see a decrease of the final particle mass, but it is bouncing 
that balances it. In the low density model, particles grow 
only in the Hit & Stick (SI) regimes. When particles leave 
these regimes, the growth stops due to bouncing. The bor- 
der of the SI regime is determined by the collision energy 
being lower than 5 x i?roii, where i?ron is the rolling en- 



ergy of monomers (see Paper I). As the collision energy is 
EcoW — 1/2/^(^1')^ J particles in strong turbulence leave the 
SI regimes at lower particle masses. 

On the other hand, the MMSN and high density mod- 
els show that the maximum mass of the particles can even 
increase with a. The precise value of the max(m) is deter- 
mined by the intensity of the peak in the mass-density plots 
fSect. H?3.ip and this may vary somewhat between the sim- 
ulations. In the 'Mtld-4mlOO' model we have argued that 
the spike is exceptionally pronounced due to the high prob- 
ability of Sticking through surface effects (S2) collisions at 
the initial part of the fast growth. However the main point 
is that in the MMSN/high density simulations the maxi- 
mum particle masses all end up around 1 g, independent of 
the turbulent strength. 

The reason for this is the nature of the S2-B2 transition, 
which occurs at projectile masses of 10"'' g in the 'pP' 
plot. As explained before, collisions in the 'pP' plot are 
the only way by which particles can grow after the Hit 
& Stick (SI) phase is finished. Thus, we require a broad 
distribution for a high growth rate. However, B2 collisions 
works in the opposite way: it transfers mass from the target 
to the projectile, narrowing the distribution and decreasing 
the overall probability for the 'pP' process. Thus, once B2 
becomes effective, there is a shift from the 'pP' panel to 
the 'pp' panel. For the MMSN/high density models this 
behavior is always present and the important quantities 
involved (i.e., relative probability of B2 over S2) scale with 
mass and not with velocity. The result is that the maximum 
masses particles achieve are ~ 1 g and rather insensitive to 
the strength of the turbulence. 

4.6. Varying the critical mass ratio 

We perform simulations in the disk models with a = 10~^ 
but with a varying critical mass ratio. We explore how the 
dust distributions change upon using r„i — 10, 100 and 
1000. Table [H lines 10 to 18, shows the parameters de- 
scribing the distribution functions, and Fig. Illb illustrates 
the maximum particle mass as a function of the critical 
mass ratio. 

By examining Table [T] we see that using = 10 in 
the low density model ('Ltld-4ml0') results in heavier and 
more compact particles. The low critical mass ratio means 
that the biggest particles in the different sized regimes can 
sweep up the projectiles and grow to bigger sizes, eventually 
reaching the fragmentation line, where growth stops. As 
discussed in Sect. 12.21 assuming a fragmentation velocity of 
1 m s"'^, the maximum Stokes number of the aggregates is 
4.7 X 10""^. This value is almost reached in this model. 

We find that there is no significant difference between 
the Tin = 100 and 1000 simulations in the low density 
model. The explanation for this can be found by exam- 
ining the width of the mass distribution in the Hit & Stick 
(SI) phase. This initial phase is happening in the same way 
independently of the critical mass ratio. If the critical mass 
ratio Vm is equal to or larger than the width of the distri- 
bution function, collisions between different size particles 
in the 'pP' regime are inhibited. After the SI phase, the 
width of the distribution in the low density regime is ap- 
proximately 100. Therefore, we do not see any difference 
when the mass threshold is shifted from = 100 to 
— 1000; in both cases collisions occur between equal-size 
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Fig. 9. Same as Fig. |3]but for the high density model. As in Fig. [5^, we zoom in on the peak at the mass distribution. 
The solid white line indicate how long our 'local approach' assumptions are valid at i = 10^ yr (discussed in Sect.|3l. 
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Table 1. Overview and results of all the simulations. 
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In this table. Col. 1 describes the model names. 'L' stands for the low density model, 'M' is the MMSN model, 'H' is the high 
density model, the letter 't' and the following number indicates the value of the turbulence parameter, the letter 'm' and the 
number shows the used critical mass ratio values. Columns 2, 3 and 4 shows the gas density, turbulence parameter and the critical 
mass ratio respectively. Columns 5, 6, 7, 8 and 9 list the parameters defined to characterize the distribution functions. These are 
the average maximum mass in Col. 5, the average maximum enlargement parameter in Col. 6, the minimum enlargement parameter 
in Col. 7, the end of the compaction phase in Col. 8 and finally the average maximum Stokes number in Col. 9. 
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Fig. 11. The maximum mean particle mass as a function of the turbulence parameter (a) and critical mass ratio (b). 



particles only and these are either SI or (when this stage is 
over) Bl. 

For the high density models (MMSN/Desch) we find 
that the outcome is again similar: growth halts at 0.1 
g (within a factor of 10) and no clear dependence on r,„ 
is seen. For the high mass ratios, growth is always in the 
similar-size regime. Here, it is the gas density that deter- 
mines the velocity, i.e., whether we have a sticking (SI) or 
a bouncing (Bl) collision. Therefore, if = 1000, the high 



density model produces heavier particles than the MMSN 
model (see Fig. Illb). For lower it is again the nature of 
the S2-B2 transition regime that limits the maximum mass. 

Thus, the critical mass ratio is an important parame- 
ter since it determines the relative likelihood of collisions 
occurring in the different-size regime, which are in general 
more conducive to growth. Conversely, in simulations where 
B2 collisions are important - which have the effect to nar- 
row the distribution ~ the width of the distribution will 
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Fig. 12. The collision frequencies of the 9 collision types in the MMSN model with a = 10^^ and 



100. 



correspond to the value of the parameter, although we 
have also seen that the absolute size/mass is rather insen- 
sitive to it. Overall, these arguments indicate that a good 
knowledge of this parameter is important. 



5. Discussion 

We performed simulations with varying turbulence param- 
eter and critical mass ratio values in three disk models hav- 
ing low, intermediate and high gas densities. We find that 
Hit & Stick (SI) and Bouncing with compaction (Bl) are 
the most dominant collision types. All simulations show the 
presence of long lived, quasi-steady states. Fragmentation 
is rarely present, but even then, only for a limited time pe- 
riod. The absence of fragmentation is due to the bouncing 
collisions. 



5.1. The sensitivity of the results 

As presented in Sect. HI the outcome of our simulations 
is determined by the collision kernel, the relative velocity 
field. A significant change in one, or both can alter the 
evolution of the aggregates. 

Here we present the results of a test simulation, where 
the Sticking through surface effects (S2) - Bouncing with 
mass transfer (B2) transition regime in the 'pP' plot is 
neglected and replaced by S2 collisions. This alternative 
transition regime provides a good opportunity to further 
examine the fast growth presented in Sect. 14.3.11 as the 
kernel is now simplified. The new kernel also gives us the 
possibility to see how much the outcome of our simulations 
can be altered by changing critical areas of the parameter 
space. As the transition regime is only constrained by one 
experiment in a rather small area (see e.g. Fig. or Fig. 11 
in Paper I), further experiments may make it necessary to 
change this part of the parameter space. We use the same 
initial conditions as in the 'Mtld-4mlOO' model described 
in Sect. 



In this case, the heaviest particle experiences increased 
relative velocities, as soon as it reaches m = 0.1 g, and 
the particle undergoes a fast growth period (as in the orig- 
inal MMSN simulation, Sect.|13]). Figure [13] illustrates the 
growth timescale of the heaviest particle (dotted line) and 
the minimum growth timescale possible (solid line). As 
there are no B2 collisions to reduce the mass of the heaviest 
particle, the growth timescale reaches the maximum that is 
possible. The heaviest particle increases in mass until the 
rest of the particle population enters the B2 regimes above 
0.1 g in the 'pP' plot. In this simulation, the maximum av- 
erage mass is 27 g, whereas in the original simulation with 
the transition regime, the value is 4.18 g. 

This work, together with Paper I, is the first attempt to 
calculate dust growth in protoplanetary disks on an empir- 
ical, thus more realistic basis. However, a few more cycles 
of the feedback loop between the laboratory experiments, 
the models of the kind described by Paper I and the models 
described in the paper have to be conducted before we can 
get near a truly reliable model of dust growth in protoplan- 
etary disks. 

5.2. Retention of small grains 



iDullemond fc DominikI (|2005l ) showed that without a 
mechanism that reduces the sticking probability of parti- 
cles in the upper layers of the disk or without a continuous 
source of small particles, the observed SEDs of TTauri 
stars would show very weak infrared exces s. The SEDs 
of TTa uri stars have strong IR excess (e.g. iFurlan et al.l 
(|2005l ). iKessler-Silacci et all (|2006l )): therefore, some kind 
of grain-retention mechanism is needed to explain these 
SEDs. Previous models of grain growth assumed a continu- 
ous cycle of growth and fragmentation, whic h provides the 
necessa ry amount of small partic l es (se e e.g. iBrauer et al.l 
(|2008ar ). IDullemond fc DominikI (|2005l ). Birnstiel et al. 
(2009)). Our simulations, however, showed that the mass 
distribution function is narrow. Small, monomer sized 
particles are not present and fragmentation is ineffective 
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Fig. 13. Growth timescale in the test simulation where 
the S2-B2 transition regime is replaced by S2 collisions 
only. The dotted line represents the growth timescale of 
the heaviest particles, the solid line is the minimum growth 
timescale. In this scenario, the growth timescale reaches the 
maximum possible value. 



in providing small particles, which could be transported to 
disk atmospheres. The question naturally arises: how can 
small grains be produced in our collision model? 

One possib l e solu tion might come from bouncing. 
IWeidling et"all (|2009l ) performed bouncing experiments 
by putting an aggregate onto an oscillating metal plate 
and measuring the porosity of particles due to collisions 
with the plate. They observed that approximately 10% 
of the projectile mass eroded during the experiment (see 
Table 1 of their paper). This mass loss can happen due to 
the initial collisions; thus the eroded mass slicked to the 
baseplate. It is also possible that small pieces of fragments 
grind off when the aggregates bounce, which cannot be 
observed in the experiment. These ground off particles can 
then diffuse out of the midplane and provide the necessary 
amount of small particles to the upper layers of the disk. 
Future laboratory experiments are needed to quantify the 
level of ground off particles in bouncing collisions. 

The second possible explanation is provided by dust 
growth at the upper layers of the disk. We performed two 
simulations at four pressure scale-heights in the low density 
model and in the MMSN model using a = 10"''. We find 
that the relative velocity of two monomers in the Brauer 
model is 2 m s~^, thus monomers at these heights do not 
coagulate, only bounce. The particles in the MMSN model 
can form aggregates of maximum of 10 /im in size. Using 
a higher a (as is mostly assumed in the upper layers of 
the disk) can completely halt even this limited growth. 
Therefore, bouncing could be the key ingredient the mech- 
anism that reduces the sticking probability of the parti- 
cles. However, if substantial vertical turbulent mixing takes 
place, this may not help, because these monomers would 
then be "vacuum cleaned" away by the bigger particles at 
the interior of the disk. Further studies of ID vertical slices 
of disk models are needed to investigate this scenario. 



5.3. Implications for planetesimal formation models 

One can also see that coagulation only cannot produce 
planetesimals with the conditions presented in this work. 
Even if the turbulence parameter is taken to be zero, 
relative velocity due to radial drift is preventing particles 
to cross the so called 'meter size barrier'. An ideal environ- 
ment for particle growth is a pressure bump in the dead 
zone where both the turbulent and radial relative velocities 
are reduced . Such an env ir onme nt is located around the 
snow line (|Kretke fc LinI (|2007l )). iBrauer et all (|2008bD 
showed that in these pressure bumps relative velocities 
stayed below a presumed fragmentation threshold of 10 
m s"'^, presenting a window through which particles can 
overcome the m-size barrier, although they assumed perfect 
sticking (no bouncing) below the fragmentation barrier. 
Future studies have to verify whether planetesimals can be 
formed with the collision model presented in this study. 

Another planetesimal forming mechani sm is the grav- 
itation al collapse of swarms of boulders ([Johansen et al.l 
(I2007D V This scenario assumes that large amount of 
the solid material is presented in dm sized boulders 
{St > 0.1) at the midplane of the disk. These boulders 
then concentrate in long-lived high pressure regions in the 
turbulent gas and these initial over-densities are further 
amplified by the streaming instability. This mechanism 
forms 100 km sized objects on a very short timescale (some 
orbits). However, our simulations produce particles with 
St « 10^** which is due to Bouncing with compaction 
(Bl) and the low (1 m s^'^) fragmentation velocity of 
silicates. Using a 'stickier' material such as ices or parti- 
cles with organic mantels may produce bigge r particles. 
Molecul ar dynamic simulati ons (e.g. Domini k fc TielensI 
([T997h . iWada et al.l (pOOTh . iWada et al.i (|200^) showed 
that icy aggregates could have fragmentation velocities 
of about 10 m s~^, although these findings have yet to 
be confirmed by laboratory experiments. Similarly, it is 
conceivable that the enhanced sticking capabilities of ices 
will prevent the bouncing, which is so omnipresent for 
small particles in our simulations, or shifts it to larger sizes. 

ICuzzi et al.l (120081 ) outlined an alternative concentra- 
tion mechanism to obtain gravitationally unstable clumps 
of particles, which can then undergo sedimentation and 
form a 'sandpile' planetesimal. In this model turbulence 
causes dense concentratio ns of aerodynamica lly size-sorted, 
chondrule-size particles (jCuzzi "eTall (l200ll )')- more pre- 
cisely, particles of Stokes numbers St = Re~^^'^ « 10~* in 
our simulations. Since growth in our models is typically 
halted at these Stokes numbers, this concentration mech- 
anism is an obvious successor to coagulation - at least 
where it concerns the conditions adopted in this paper 
(lAU, sihcates). 

However, it should be emphasized that the formation 
of a gravitationally unstable clump does not imply plan- 
etesimals will form unimpededly. An important question 
to address is h o w co llisions will affect the collapse. In 
the[C uzzi "eTall (pOOSh scenario the collapse occurs on a 
sedimentation timescale and for these high densities col- 
lisions between particles will be frequent. Likewise, in the 
Johansen scenario - where the collapse occurs on an orbital 
timescale and involves St ^ 0.1 particles - collisions can 



22 



A. Zsom et al.: The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals? II. 



be rather violent. CoUisional fragmentation or erosion may 
change the appearance of the collapse, because the small 
fragments are carried away by the gas. The role of collisions 
in these situations is certainly an important question, and 
our new collision model provides a tool to quantitatively 
address this issue in future studies. 

5.4. Consequences for laboratory experiments 

One can see from Fig. 11, in Paper I, that only a small part 
of the parameter space is covered by experiments. Although 
laboratory experiments cannot be made at every point of 
the parameter space, we suggest future ones based on Figs. 
[5l [7] and [TU] in order to better understand dust growth in 
the early stages of planet formation. 

— More experiments in the 'cc' and 'cC regimes are 
needed as particles get compactified by the end of their 
evolution. Thus, most of the collisions happen in this 
regime, at velocities between 0.1 and 100 cm s~^, at 
masses between 10~^ and 10 g. 

— As seen in Figs . [5l [7l and fTOl the 'hot spots', where most 
of the collisions are happening, are located in the equal 
sized regimes, at the left side of the fragmentation line. 
Therefore, it is important to map these areas of the 
parameter space in detail. 

— We define a sharp border line between the Hit & Stick 
(SI) and Bouncing with compaction (Bl) collisions. If 
there is a continuous transition between SI and Bl, the 
growth of particles would not be halted by bouncing at 
such low particle sizes. As many collisions are happening 
in the 'pp' and 'cc' regimes, even a small probability of 
growth could increase the particle sizes. 

— As seen in Fig.|6b, particles in high gas density environ- 
ments can have enlargement parameters much higher 
than 6.6 (0 = 0.15). An interesting question is whether 
the collision types and regimes are also valid for parti- 
cles with such a low volume filling factors, or whether 
these particles have a different collision behavior? 

— The Sticking through surface effects (S2) - Bouncing 
with mass transfer (B2) transition regime greatly affects 
the outcome of the simulations (see Sect. 15. ip . However, 
the transition regime is only mapped at the high velocity 
and low mass regions. Therefore, it is essential to better 
constrain this part of the parameter space. 

— The critical mass ratio affects the particle masses and 
porosities. Experiments are needed to constrain its 
value. 

— The bouncing model, described in Paper I, has impor- 
tant implications for the evolution of dust aggregates in 
protoplanetary disks but it is unfortunately still based 
on too few experiments. Further experiments are needed 
to refine the model, as Bouncing with compaction (Bl) 
is the most frequent collision type in all of the simula- 
tions. 

6. Summary 

We performed simulations of dust growth using the Monte 
Carlo code of ZsD08 and a dust collision model based on 
laboratory experiments (Paper I). We performed simula- 
tions at the midplane of three disk models having low 
(2.4 X 10"-^^ g cm^'^), intermediate (1.4 x 10~^ g cm~^) and 
high (2.7x 10~^ g cm~^) gas densities at 1 AU distance from 



the central star. We vary the turbulence parameter (a) and 
the critical mass ratio (r^) to explore their effects on the 
mass and porosity distribution functions. Our main results 
are: 

— Upon using a = 10"", the low density / MMSN / high 
density model produces particles with maximum mean 
mass of 9.7 x 10~^ g / 4.18 g / 0.23 g, the maximum av- 
erage enlargement parameter of these particles are 7.12 
/ 21.9 / 38.0. The maximum average Stokes numbers 
are 2.2 x lO^"* / 2.8 x 10^* / 3.6 x 10"^. 

— We find that particle evolution does not follow the pre- 
viously assumed growth-fragmentation cycles. Although 
catastrophic fragmentation is present for a short pe- 
riod of time in some of the models (typically when 
a = 10^'^), it has a fringe effect. Particles in most of 
the simulations do not reach the fragmentation barrier 
because their growth is halted by bouncing. 

— We see long lived, quasi-steady states in the distribution 
function of the aggregates due to bouncing. The final 
equilibrium state is not reached within the drift or the 
viscous timescales. 

— We performed simulations with varying turbulence 
strength. We find that the system is 'non-linear': The 
maximum mass of particles is not a decreasing function 
of the turbulence parameter and is not an increasing 
function of the gas density. 

— We explored the effects of the critical mass ratio. We 
find that different critical mass ratios can affect the 
particle evolution. Small critical mass ratios can pro- 
duce heavier particles, while big values of can halt 
the growth earlier. 

— The maximum Stokes number is rather independent of 
the gas density and the strength of the turbulence. 

— The maximum mass of the aggregates is limited to 1 
g due to the S2-B2 transition regime. 

— The Stokes number lO"** particles can be concentrated 
in turbulence by aerodynamical size-sorting, thus plan- 
etesimals can form from these particles. 
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Appendix A: Time evolution and animations 

We offer six animations to illustrate the time evolution 
of the aggregates. These animations can be viewed at 
www . mpia . de/homes/zsom/Site/ SLnimat ions2 . html. The 
mass_vs_psi_lowd.mpg file shows the time evolution of 
the masses and enlargement parameters of our repre- 
sentative particles. The x-axis is the particle mass in 
gram units, the y-axis is the enlargement parameter. The 
regimes-lowd.mpg file illustrates the time evolution of 
Fig. [S] The contour levels are collision frequencies / grid 
cell. The same movies are provided for the MMSN model, 
these are mass_vs_psi_mmsn.mpg and regimes_mmsn.mpg 
respectively. And also for the high density model: 
mass_vs_psi_highd.inpg and regimes Jiighd.mpg. 
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